UN  I 


HARVARD    UNIVERSITY 

Library  of  the 

Museum  of 

Comparative  Zoology 


OCCASIONAL  PAPERS  ,,,^ 

of  the  ' 

MUSEUM  OF  NATURAL  HISTORY 

The  University  of  Kansas 
Lawrence,  Kansas 


''^'lo 


NUMBER  108,  PAGES  1-9  29  FEBRUARY  1984 


ESTIMATION  OF  THE  LATITUDE  AND  LONGITUDE 
COORDINATES  OF  POINTS  ON  MAPS 

By 

C.  HaMAKER'    and  J.  W.  KOEPPL^ 

Considerable  information  on  the  distribution  of  plants  and  animals  is 
represented  as  points  on  maps.  For  example,  presence  of  a  particular 
taxon  is  usually  indicated  by  a  dot  or  other  symbol  at  the  collection 
location  (Fig.  1;  Duellman,  1970;  Hall,  1981;  Hulten,  1950;  Lee  et  al., 
1980;  many  others).  Plotting  many  such  occurrences  on  a  map  yields  an 
estimate  of  the  geographical  range  of  a  taxon.  Placement  of  the  points  is 
generally  determined  from  written  descriptions  of  locations  provided  by 
the  collector.  Plotting  many  such  locations  is  laborious  and  time  consum- 
ing. When  comparison  or  consolidation  of  the  distributions  on  several 
maps  is  desired,  the  process  is  hindered  by  the  different  map  scales  and 
projections  where  this  information  is  found  and  by  the  large  number  of 
points.  This  process  can  be  expedited  by  rendering  all  point  locations  into 
a  common  spherical,  latitude  and  longitude  coordinate  system  and  em- 
ploying computerized  cartographic  methods  to  display  the  distributions  on 
a  new  map  of  appropriate  scale  and  projection.  Herein,  we  describe  two 
general  methods  for  translation  of  point  locations  into  latitude  and 
longitude  coordinates.  We  test  the  methods  on  several  commonly  em- 
ployed map  projections  and  discuss  the  characteristics,  strengths  and 
weaknesses  of  the  methods.  Moreover,  because  maps  are  merely  sets  of 
points,  and  the  methods  described  here  are  general,  the  methods  have 
potential  for  translating  whole  maps  from  one  projection  to  another. 

METHODS 

The  basic  situation  presents  a  point  P,,  with  known  Cartesian  coordi- 
nates (Xq,  yo)  and  unknown  spherical  coordinates  {slq,  by)  along  with  a 


'  Department  of  Mathematics,  The  University  of  Kansas,  Lawrence,  Kansas  66045. 
-  Museum  of  Natural  History.  The  University  of  Kansas,  Lawrence,  Kansas  66045. 


OCCASIONAL  PAPERS  MUSEUM  OF  NATURAL  HISTORY 


5  5'N  6  0*.\ 


"^zr 


,JV_ 


""?■     /  ARCTIC  OCEAN  ^  X'i-  ^ 

>     1^     CHUKCHI  SEA  ^    "J-X  i 

BERING  SEA       '\       T^  (  .  ■•  S  V:?^  <(    ^v^Q^i'^-^ 


3* 


PACIFIC  OCEAN 


NORTHWESTERN    NORTH    AMERICA 
LAMBERT   CONIC   PROJECTION 
SCALE    I  :  18.000,000 


J 

!o  J  HUDSON 
^  ^      BAY       ^ 


a,.-.=Q^i.  J., 


Figure  1  .—Example  of  a  dot  distribution  map  of  Spermophilus  parryii  taken  from  Pearson, 
1981.  Known  spherical  coordinates  may  be  obtained  from  intersection  of  latitude  and 
longitude  lines,  distinctive  geographic  or  political  boundaries  or  locations  of  towns  and  cities. 
Corresponding  Cartesian  coordinates  obtained  from  a  digitizer  provide  the  other  set  of 
coordinates  comprising  the  reference  coordinate  system.  Cartesian  coordinates  of  dots  are 
then  translated  into  their  spherical  coordinates  using  the  methods  described  in  this  paper. 


collection  of  points  P^,,  k=  1 n,  whose  Cartesian  coordinates  (x^., 

y^_)  and  spherical  coordinates  (a^^,  b,^)  are  both  known.  The  problem  is  to 
estimate  (a,,,  bo).  We  assume  that  the  transformation  (inverse  map 
projection) 


(a,b)  =  [a(x,y),  b(x,y)]  =  F(x,y) 


(1) 


from  Cartesian  to  spherical  coordinates  is  smooth.  Because  the  estimate 
of  (ay,  by)  is  to  be  made  without  knowing  the  transformation  F.  we 
devised  two  types  of  approximations  to  F  by  least  squares  techniques. 
Global  method  of  order  r.— Here  (a,b)  =  F(x,y)  is  approximated  through 

the  region  containing  P, P„  by  a  polynomial  of  specified  degree  r, 

i.e. 


a=     E 


m 


q: 


and 


m  =  0  j  =  0 

r        m 

b=    E      s:   iS 

m  =  0  j  =  0 


mj 


m-j 


inj 


xJ  y 


XJ  ym-J 


(2) 


(3) 


The  coefficients  q:,„|  and  /S,,,,  are  found  by  standard  least  squares  tech- 
niques, described  below  (see  Strang,  1980). 

Define  the  matrix  M  and  vectors  W,  Z.  A,  and  B  by 


ESTIMATION  OF  POINTS  ON  MAPS 


1=  1  +- +  J 

2 


Mk,,  =  (Xk>i(yk)'^-J,  where 


0  ^  m  ^  r  (4) 


1  ^  k  ^  m 

W.  =  a„,j,       Z,  =  i3^j  (5,6) 

Ak^a^,  and  Bk  =  bk  (7,8) 

In  this  notation,  substitution  of  the  respective  coordinates  of  P,,  .  .  .,  ?„ 
into  (2)  and  (3)  gives  the  matrix  expressions  A=^MW  and  B  =  MZ.  The 
best  choices  of  W  and  Z  in  the  sense  of  least  squares  satisfy  the  normal 
equations 

MTMW  =  MTA  where       M^  is  the  transpose  of  M;  (9) 

MTMZ  =  MTB.  (10) 

After  solving  (9)  and  (10)  for  the  coefficient  vectors  W  and  Z,  an 
approximation  of  (slq,  bo)  can  be  computed  by  using  (x,y)  =  (Xo,yo)  in  (2) 
and  (3). 

Local  method  of  order  r.— The  global  method  described  above  approxi- 
mates the  transformation  F  thoughout  the  map  by  using  all  known  points. 
In  contrast,  the  local  method  approximates  only  (ao,bo)  =  F(Xo,yo)  while 
weighting  the  known  points  P^.  closest  to  Pq. 

By     translation     of    Cartesian     coordinates,     we     may     assume 
(Xo,yQ)  =  (0,0).  Taylor's  Theorem  states  that 

Fi(x,y)  =  a(x,y)=     E        E    a^j  xj  y^-J +  E(x,y)  (11) 

m=Oj=0 

where  the  polynomial  of  degree  r  approximating  F,  has  coefficients 
Note  that  ao  =  a(0,0)  =  ao,o-  The  error  term  satisfies 

r+l 

|E(x,y)l^Cr(x2  +  y2)  2  (13) 

where  the  positive  constant  C^  depends  on  the  maximum  values  of  the 
(r+  l)-order  partial  derivatives  of  F,  in  the  region  of  interest.  If  these 
partials  are  small,  then  multiplying  (11)  by 

r+l 
d,  =  (x2  +  y2)      2  (14) 

gives  an  approximation  with  an  error  bounded  independently  of  (x,y)  by 
the  small  constant  C^: 

r        m 
dra(x,y)=     E        E    ttmj  dr  xJ  ym-j     .  (15) 

m=Oj=0 
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Using  the  coordinates  {x^,,y^.)  and  a^^,  of  the  reference  points  in  (15) 
yields  the  linear  system: 

MZ^A  (16) 

where  the  matrix  M  and  the  vectors  Z  and  A  are  defined  as  follows.  Let 


drk  =  (^k-  +  yr) 


r+l 

2 


and 


Mkj  =  drk  (Xk)J(yk)'"-J  where 
Ak  =  drkak. 


i=l  + 


m2  +  m 


+  J 


O^m^r 
O^j^m 
0<k<n 


(17) 


(18) 


(19) 


(20) 


As  in  the  global  method,  the  best  choice  of  Z  in  the  sense  of  least 
squares  is  found  by  solving  the  normal  equation.  However,  the  only  entry 
in  Z  which  we  require  is  Z,  =aQQ  =  ao,  which  can  be  found  by  Cramer's 
Rule  or  by  a  partial  Gaussian  elimination.  The  final  formula  via  Cramer's 
Rule  for  the  approximation  of  ag  is 


ao  = 


det(R.,0 
detCSuO 


(21) 


where,  for 


1=  1  +^-^  +  m. 


O^j^r 
O^m^j 


(22) 


and 


0     i,P"  +  P,  u  0<p<r 

F=l+      -^+q,  where  r^  ~ 

2         ^  0<a<n 


o^q^p' 


(23) 


we  set 


n 


S,,=    E     d,k^(Xk)i  +  q(yk)-  +  P-J-i 
k=l 


(24) 


and 


R,f  = 


L     ak  drk  (Xk)q(yk)P-H     iff=l 
k-1 


S,r 


if  f>l 


(25) 
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The  approximation  of  b,^  is  obtained  by  applying  the  same  considera- 
tions to  the  Taylor  expansion  of  F-,(x,y)  =  b(x.y).  The  final  formula  is  the 
same  except  that  b^^  is  substituted  for  a,^  in  the  definition  of  R,,'. 

IMPLEMENTATION  AND  SIMULATION 

The  methods  described  above  are  routinely  implementable  in  a 
FORTRAN  program.  We  supplemented  that  program  with  other  sub- 
routines in  order  to  simulate  the  basic  problem  and  test  the  effectiveness 
of  these  methods. 

1)  Spherical  coordinates  for  the  set  of  reference  points  in  the  interval 
1-89  degrees  latitude  and  1-179  degrees  longitude  were  supplied  by  a 
random  uniform  number  generating  subroutine.  The  set  of  nine  test  points 
were  uniformly  distributed,  having  spherical  coordinates  (15°. 30°), 
(15°,90°),  (15°, 150°),  (45°, 30°),  (45°,90°),  (45°, 150°),  (75°, 30°), 
(75°,90°),  and  (75°, 150°). 

2)  For  each  of  four  map  projections,  Cartesian  coordinates  were 
computed.  Our  subroutines  used  the  equations  provided  by  Maling  (1973) 
for  Lambert  Azimuthal  equal-area  projection,  Mercator  projection,  Si- 
nusoidal (Sanson-Flamsteed)  projection,  and  Hammer-Aitoff  projection 
(Fig.  2). 

3)  For  each  test  point,  its  Cartesian  coordinates  and  the  Cartesian  and 
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Figure  2.  Four  projections  for  which  simulations  were  conducted:  a,  Mercator  projection;  b, 
Azimathal  equal-area  (Lambert)  projection;  c.  Sinusoidal  (Sanson  Flamsteed)  projection;  d, 
Hammer-Aitoff  projection.  Simulations  were  performed  only  in  the  first  quadrant.  The  nine 
points  on  each  projection  were  the  known  points  whose  spherical  coordinates  were  estimated. 
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spherical  coordinates  of  its  nearest  neighbors  among  the  reference  points 
were  used  to  approximate  its  spherical  coordinates.  The  number  of 
nearest  neighbors  used  was  varied  for  each  point,  as  were  the  order  r  of 
the  method  and  the  total  number  of  reference  points.  These  computations 
were  performed  for  each  map  projection  and  for  both  the  local  and  global 
methods. 

4)  As  a  rough  measure  of  the  effectiveness  of  these  methods,  we 
computed  the  variance  of  the  approximated  spherical  coordinates  from 
the  actual  spherical  coordinates  of  the  test  points. 

RESULTS 

Results  of  the  simulations,  when  varying  the  order,  number  of  nearest 
neighbors,  density  of  the  reference  coordinates  and  method  of  computa- 
tion, are  reflected  in  the  variances  (Table  1).  The  variance  (s-)  was 
computed  as  the  sum  of  the  squared  deviations  of  latitude  and  longitude 
(djTj  +  djT^)  divided  by  the  number  of  determinations  (n): 

n 

S2=    E      (d2^  +  d2g)/n 
i  =  1 

A  variance  of  less  than  1  is  attained  when  the  deviations  of  the  latitude 
and  longitude  are  both  less  than  2  "  '^'  or  one  is  0  and  the  other  is  less  than 
1. 

For  all  the  simulations,  increasing  the  order  of  the  computations 
markedly  decreased  the  variance  and  hence  the  error  between  the  actual 
and  estimated  spherical  coordinates;  increasing  the  number  of  nearest 
neighbors  for  any  order  increased  the  variance  (Table  1). 

The  number  of  random  reference  coordinates,  uniformly,  but  ran- 
domly distributed  in  the  first  quadrant  of  the  map  (Fig.  2)  was  either  78 
or  275.  These  numbers  may  be  regarded  as  densities,  that  is,  the  number 
of  points  per  unit  area.  Together  with  the  number  of  nearest  neighbors, 
these  two  values  are  reflected  in  the  mean  distance  between  the  point 
whose  spherical  coordinates  are  to  be  estimated  and  the  reference 
coordinates  used  in  the  estimation  (Table  1).  Comparison  of  variances  of 
the  four  projections  show  that  results  vary  depending  on  the  projection 
and  the  parameters  chosen  (Table  1).  No  single  projection  is  best  over  all 
ranges  of  parameters.  A  similar  result  emerges  when  we  consider  the 
global  versus  the  local  methods.  The  exception  to  this  finding  is  that  the 
local  method  is  more  robust  when  all  available  coordinates  are  used  to 
form  the  estimate. 

DISCUSSION 

For  either  the  global  or  the  local  method,  if  our  sole  goal  is  to 
estimate  the  spherical  coordinates  to  a  high  accuracy,  and  other  consid- 
erations are  only  minor,  using  high  order  fits  will  achieve  satisfactory 
results.  In  theory,  there  is  no  limit  to  the  accuracy  that  can  be  achieved. 
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The  only  difference  between  the  two  methods  is  that  they  will  achieve  the 
desired  accuracy  at  somewhat  different  rates  depending  on  a  variety  of 
factors  already  mentioned,  and  some,  such  as  error  in  the  coordinates 
supplied,  which  we  will  consider  briefly  later. 

Unfortunately,  there  are  a  number  of  practical  limitations  which  bear 
on  the  actual  accuracy  that  can  be  achieved. 

The  order  of  the  computations  possible  is  limited  by  the  precision  of 
the  computer  and  the  magnitude  of  the  numbers  that  can  be  processed. 
The  Honeywell  66/60  on  which  the  simulations  were  performed  allowed 
computations  in  the  range  ±  10-''*^.  We  found  that  we  could  compute  5th 
order  global  and  3rd  order  local  estimates  without  encountering  size 
errors. 

The  number  and  distribution  of  the  reference  coordinates  can  greatly 
influence  the  results.  The  more  reference  coordinates  distributed  through 
the  area  in  which  the  points  are  to  be  estimated,  the  better  the  estimates 
tended  to  be.  Most  of  these  coordinates  must  be  determined  and  input 
prior  to  the  analysis;  this  may  be  very  time-consuming.  It  is  better  to  limit 
the  number  of  reference  coordinates  for  this  reason.  On  the  other  hand, 
we  have  seen  that  the  miminum  number  of  coordinates  for  any  order  (r) 
of  determination  is  limited  by  the  expression  (r+ l)(r  +  2)/2.  Hence 
algorithms  of  orders  1,2,  and  3,  require  a  minimum  of  3,  6,  and  10 
reference  coordinates,  respectively.  In  practice,  it  is  best  to  use  more 
reference  coordinates  than  this  because  certain  distributions  of  reference 
coordinates  produce  ill-conditioned  matrices  in  the  normal  equations.  In 
theory,  such  a  matrix  need  not  be  invertible  (e.g.,  if  the  known  points  P^, 
all  have  the  same  latitude).  Thus,  we  must  assume  that  the  points  P^^  are  in 
"general  position,"  and,  as  a  rule  of  thumb,  we  use  (r-l-2)(r-l-3)/2 
reference  coordinates,  the  minimum  number  theoretically  required  for 
order  r  +  1 ,  when  computing  a  fit  of  order  r. 

The  number  of  nearest  neighbors  used,  besides  agreeing  with  the 
constraints  of  minimum  sample  size  and  number  of  reference  coordinates, 
depends  on  the  type  of  projection  and  map  scale.  Large  scale  maps  which 
use  projections  such  as  the  sinusoidal,  having  singularities  or  areas  of 
severe  distortion  in  the  area  of  interest,  are  best  limited  to  fewer  nearest 
neighbors,  while  small  scale  maps  of  projections  having  few  singularities 
such  as  the  Mercator  and  more  gradual  distortions  may  use  larger 
numbers  of  nearest  neighbors. 

A  similar  choice  exists  when  considering  global  versus  local  al- 
gorithms. Global  fits  are  best  used  for  map  projections  and  scales  with 
less  severe  distortion,  while  the  local  fit  is  better  for  projections  and 
scales  where  more  distortion  is  present.  The  choice  of  number  of  nearest 
neighbors  and  local  or  global  algorithms,  to  some  extent,  mediates  the 
practical  limits  of  computing  fits  of  high  order. 

Our  simulations  tested  the  methods  discussed  while  using  error  free 
data.  We  did  this  by  computing  the  Cartesian  coordinates  to  high 
precision  for  spherical  coordinates  already  known.  In  practice,  the 
researcher  will  introduce  additional  error  due  to  imperfect  measurement 
of  coordinates.  We  have  not  studied  how  our  computations  are  affected  by 
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such  errors  in  measurement.  We  feel  that  the  error  of  point  determination 
by  the  local  method  will  prove  to  be  proportional  to  its  weighting  factor, 
while  the  global  methods  error  will  be  distributed  more  evenly  over  the 
reference  coordinates  considered. 

In  the  final  analysis,  one  faces  a  series  of  compromises  or  trade-offs 
with  respect  to  the  parameters  chosen  for  the  estimation  of  spherical 
coordinates.  The  choices  one  makes  will  be  determined  by  the  problem  at 
hand,  and  the  resources  and  requirements  of  the  researcher.  In  our  view 
none  of  these  considerations  severely  limit  the  utility  of  the  methods 
presented. 
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SUMMARY 

To  estimate  the  latitude  and  longitude  coordinates  or  arbitrary  points 
on  maps,  we  select  reference  points  whose  coordinates  are  known  in 
Cartesian  and  spherical,  latitude  and  longitude  coordinate  systems.  Func- 
tions relating  these  two  sets  of  reference  points  may  be  approximated  and 
used  to  translate  additional  points  from  Cartesian  to  latitude  and  longitude 
coordinates.  The  methods  may  also  be  used  to  translate  maps  from  one 
projection  to  another. 
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